Temperature dependence of normal mode reconstructions of protein dynamics 
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Normal mode analysis is a widely used technique for reconstructing conformational changes of 
proteins from the knowledge of native structures. In this Letter, we investigate to what extent 
normal modes capture the salient features of the dynamics over a range of temperatures from close 
to T = to above unfolding. We show that on the one hand, the use of normal modes at physiological 
temperatures is justified provided proteins are cooperative. On the other hand, it is imperative to 
consider several modes in order to eliminate the unpredictable temperature dependence of single- 
mode contributions to the protein fluctuations. 

PACS numbers: 87.15.-v, 87.10.Tf, 87.14.E- 



The effectiveness of proteins at performing their func- 
tions strongly depends on their ability to dynamically 
explore different conformations. In many cases, the dif- 
ferent functionally relevant conformations correspond to 
global deformations of the molecules, although not nec- 
essarily of large amplitude. In such circumstances, the 
small-amplitude dynamics can be adequately described 
through Normal Modes (NM) analysis in order to deci- 
pher the structure-dynamics- function relation P, Q . 

NM analysis of macromolecules has a venerable history, 
dating back to the early 1980's and has enjoyed sev- 

eral renaissances, most prominently in its latest variant, 
the Anisotropic Network Model and its derivatives [^, @] . 
The cornerstone of NM approach is the assumption that 
functional dynamics of proteins can be captured by the 
lowest-frequency harmonic normal modes as long as the 
native-state deformations are of small amplitude. There- 
fore, the complex potential energy ruling protein 
dynamics can be expanded, as a function of the atomic 
coordinates r^, i = 1, . . . , iV, to second order in the fluc- 
tuations around the native state 

V{{r}) ^ Ar,„Ar,^, (1) 

id 

where Ar^Q, = — Ri^ , Ria being the native state coor- 
dinates. The dynamical matrix of the system is defined as 
D = M-^/'^HM-^/'^, with M the diagonal mass matrix 
and H the Hessian of the potential. The eigenvectors i/'fc , 
/c = 1, 3A^, of D are the normal modes of the protein. 

Although NM analysis has been applied with success 
also to predict large scale deformations, its appropriate 
context remains the small-amplitude regime and extrapo- 
lation to larger fluctuations should be rigorously checked 
to avoid arbitrary interpretations. In fact, at T = 300 K, 
a mere 10 % below typical protein unfolding tempera- 



tures [Tf), displacements well beyond the harmonic ap- 
proximation must be expected. A possible way to find 
suitable collective coordinates to describe large fluctua- 
tions is to use the "essential modes", the eigenvectors 
of the covariance matrix of the displacements from the 
native state, computed in a molecular dynamics simula- 
tion [7] . Unfortunately, the convergence of the covariance 
matrix toward a stationary solution is very slow Q. 

In this Letter, we establish on a firm basis the valid- 
ity of the NM description of protein dynamics and unveil 
some remarkable subtle features of this approach as tem- 
perature is increased from very low to beyond Tf. 

In order to measure the spectral weight of protein fluc- 




FIG. 1: Cartoon representations of three of the analyzed 
structures, along with the corresponding PDB codes. 
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tuations as a function of temperature, we introduce the 
thermal involvement coefhcients (TICs) 0fe(T) [9]. They 
are defined as efc(r) = {Yk{T))j^, with 



dihedral potential energy 



Yk{T) 



{{MT)-^kf)t 



(2) 



where {A(r)} = {T*(T)} - {R} is the 3iV-dimensional 
deviation of the protein structure from the native fold 
configuration projected onto the NMs, and averaged both 
over time (•••)* and over M different initial conditions 
drawn from the same equilibrium distribution at tem- 
perature T, (. . . We have exphcitly verified that the 
two averages commute, so as to strengthen the statis- 
tical significance of the sampling and averaging proce- 
dures. Statistical uncertainties were computed as the 
standard errors on the realization averages A0fc(T) = 

Protein trajectories are calculated within the isoki- 
netic scheme [lOj, which provides a correct sampling of 
the configuration space [19|. Forces are calculated using 
the coarse-grained Go model introduced in Ref. [ill where 
each amino-acid is replaced by a bead with the average 
amino-acid mass, whose equilibrium position coincides 
with that of its a-carbon. Successive beads along the 
chain are connected by stiff harmonic springs, mimicking 
the peptide bond and maintaining the chain connectivity. 
In line with native-centric schemes, non-bonded interac- 
tions between non consecutive a-carbons are modelled 
with Lennard- Jones 12-10 potentials if the atoms are in 
contact in the native state according to a given interac- 
tion cutoff Rc and with purely repulsive interactions oth- 
erwise. The parameters of the non-bonded interactions 
are fixed as Rc = 6.5 A and Rc = 7.5 A for N < 100 and 
N > 100 respectively. The force field is completed by 
standard harmonic angle-bending interactions, plus the 



N-2 
i=3 



Vdh = > ; fc^ll - cos(A0,)] + kf\l " cos3(A</.,)] (3) 



where A(j)i are the the dihedral angle deviations from the 
native values [2^ . Normal modes are computed from the 
dynamical matrix of the Go model, according to equa- 
tion ID), after a conjugate-gradient relaxation of the PDB 
structures, necessary so to obtain the true minimum of 
the model potential [2l|. 

In order to locate the unfolding temperatures, we have 
performed thermal unfolding simulations coupled to mul- 
tiple histogram reconstructions [I^ of the specific heat. 

We have analyzed a number of protein structures [23 | , 
ranging in length from TV = 57 to = 171 and het- 
erogeneous for a and P content. The first notable find- 
ing is that a restricted number of low-frequency normal 
modes is always able to capture a sizable fraction of 
the equilibrium fluctuations over a broad temperature 
range, from T = to about 10-20 % below T/. This 
is clearly represented in Figl5] showing the behavior of 
the first ten TICs of Lysozyme and of OMPR DNA- 
binding domain. The TICs of Lysozyme, after a major 
rearrangement close to T = 0, appear strikingly robust 
over the whole range where the protein structure is sta- 
ble, while those of the DNA-binding domain fluctuate 
wildly. Remarkably, in the latter structure, the mode 
whose TIC is the most important at intermediate temper- 
atures gave nearly no contribution to fluctuations close to 
the native minimum. Conversely, the NMs that captured 
the T « fluctuations appear to swap and rapidly lose 
spectral weight as temperature increases. Note that the 
first TICs of Lysozyme fall off much closer to Tf than 
those of the DNA-binding domain (see values of Tf in 
the caption). These results manifestly warn on the use 



3 




FIG. 3: Cumulative involvement coefficient over the first three modes for: (a) Lysozyme, A'^ = 162, PDB code 166L (empty 
circles) and a surrogate self-avoiding polymer occupying the same volume (filled circles) and (b) the peptide-binding fragment 
from Hsp40 (b), A'' = 171, PDB code HDJl. Insets: specific heat (unfolding) on the same temperature scale. In panel (a): wild 
type molecule (dashed line) and the surrogate (solid line). 



of single-mode reconstructions of conformational changes 
using zero-temperature NMs, as these are not guaran- 
teed to provide stable spectral measures at physiological 
temperatures. Indeed the individual low-frequency TICs 
are observed to fluctuate substantially, or even exchange 
spectral weight among each others, with respect to zero 
temperature before stabilizing. 

Despite the observed fluctuations of TIC values at low 
and intermediate temperatures, a global measure gath- 
ering the contributions from several low-frequency NMs 
appears to be a more stable indicator of spectral weight 
over the whole thermal span (see Figs. [2] (b) and [3|) 
for all analyzed structures. Note that this is also the 
case for the wildly fluctuating TICs of the OMPR DNA- 
binding domain (see again Fig. [2] (b)), when the first 
three of them are aggregated. Hence, subspace-based, 
rather than single-mode-based reconstructions carry the 
relevant spectral weight. Incidentally, this sheds further 
light on the observed correlation between low-frequency 
subspaces of NMs and essential modes pointed out in 
Ref . M- 

As Tf is approached, a transition of variable width to 
equipartition is invariably observed. Above Tf, TICs are 
no longer sensitive to structural rearrangements. How- 
ever, the nature of such transitions appears rooted in the 
special structural features of proteins. More precisely, 
TICs may relax to equipartition more or less abruptly. 
In this regard, crucial issues seem to be (i) the degree 
of structural arrangement at the secondary level and (ii) 
deviations from globularity of the scaffold. 

(i) To address the first point, we generated for a few 
structures self-avoiding polymers of the same length as 
the corresponding proteins and confined in the same en- 
veloping ellipsoid as calculated through the correspond- 



ing inertial tensor. The necessary excluded volume con- 
straints and local persistence length properties were im- 
posed by requiring that the three-body radii be never 
lower than 2.7 A, as described in Ref. 13 By doing 
so, the ensuing structures show the same internal static 
two- and three-body correlations as the original proteins, 
while almost lacking the distinctive occurrence of a and 
(3 motifs. As shown in Fig. [3l^a), NMs of such surrogate 
structures do not capture the dynamical fluctuations to 
the same extent as in the original proteins. This has been 
verified for different proteins in our ensemble. Hence, we 
reaffirm the principle that functional motions that can be 
captured through NM analysis are rooted in the special 
arrangements of protein folds ^, 15, 1^, 17 1 . This is the 
second important point we wish to stress. 

(u) Not all proteins possess globular folds, some DNA- 
and peptide- binding fragments, for example, being char- 
acterized by less regular shapes, necessary to adapt to 
their targets upon binding. This, in turn, implies a 
varying degree of structure-related cooperativity in the 
dynamics. Therefore, it is instructive to assess the im- 
pact of deviations from globularity on the temperature 
trend of TICs in such cases. In Fig. [3l^b), we show that 
the TICs describing fiuctuations of the human Hsp40 
peptide-binding fragment, a strongly non-globular pro- 
tein (see cartoon in Fig. [T]), are indeed characterized by 
a weak degree of cooperativity. In particular, it is clear 
that they no longer represent a faithful measure of the 
overall dynamics already well below Tf . 

Finally, to highlight the role of the energy landscape 
features on the temperature dependence of TICs, we have 
analyzed Fl-ATPase, a globular protein with a high con- 
tent of of a and p motifs. Fig. d] shows the TIC ther- 
mal plot computed both with the full potential and by 
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a crucial liaison between the meaningfulness of a NM- 
based description of functional fluctuations at physiolog- 
ical temperatures and the degree of cooperativity of pro- 
teins. 

The authors wish to thank Bingdong Sha for making 
the HDJl structure available prior to publication. 
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FIG. 4: Cumulative involvement coefficient over the first three 
modes for the e-subunit of the Fl-ATPase {N = 138, PDB 
code IBSN) computed with the full potential (filled circles) 
and without dihedral forces (empty circles). Lines are just 
guides to the eye. Inset: specific heat (unfolding) for the full 
model (dashed line) and without dihedral terms (solid line) 
on the same temperature scale. 



switching off the dihedral terms ([3]), which are essential 
in order to reproduce the cooperativity of protein folding 
thermodynamics within a Go scheme |l8| . Since it is also 
known that less globular proteins are less cooperative, we 
expect that, due to the lack of the dihedral term, TICs 
will also be affected, showing a behavior similar to what 
observed for non-globular proteins such as HDJl. Indeed, 
Fig. m documents the validity of our inference. Note that 
the strong cooperativity and thermal resilience of TICs 
are both profoundly eroded starting from relatively small 
temperatures. This completes the series of our findings. 

In this Letter we have shown that sizable values and 
thermal resilience of TICs are deeply rooted in the pres- 
ence of secondary motifs and in the degree of their pack- 
ing within protein scaffolds. The fluctuations of more 
irregular, less compact and less cooperative proteins are 
manifestly harder to capture through normal modes. Al- 
though our results confirm the validity of NM analysis 
of protein functional motions, we have shown that sub- 
spaces, typically spanned by the first few low-frequency 
modes, rather than single normal modes, ought to be 
considered as identifying the relevant functional mo- 
tions. This conclusion stems in a clear fashion from 
the observation that the spectral capacity of individual 
modes, as measured through TICs, often differs sub- 
stantially at physiological temperatures with respect to 
the traditional static involvement coefficients at zero- 
temperature. Moreover, we have brought to the fore 
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Since the isokinetic scheme conserves both linear and an- 
gular momenta, no special alignment procedure is needed 
before computing conformational changes, if care is taken 
to set both total momenta to zero at the beginning of 
each simulation run. 

The dihedral angle at site i is defined by the two adjacent 
planes formed by the four consecutive Ca 's at i — 2, i — 1, 

1, i + 

Codes available from the authors upon request. 

PDB codes: 166L, IFAS, lOOI, HDJl, IBFE, IBSN, 

ICHN, IFVQ, lOPC, lAWO, ITIT, lUBI, INEB. 



